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Abstract 

Data assimilation is the task to combine evolution models and observational data in order to produce 
reliable predictions. In this paper, we focus on ensemble-based recursive data assimilation problems. Our 
main contribution is a hybrid filter that allows one to adaptively bridge between ensemble Kalman and 
particle filters. While ensemble Kalman filters are robust and applicable to strongly nonlinear systems even 
with small and moderate ensemble sizes, particle filters are asymptotically consistent in the large ensemble 
size limit. We demonstrate numerically that our hybrid approach can improve the performance of both 
Kalman and particle filters at moderate ensemble sizes. We also show how to implement the concept of 
localization into a hybrid filter, which is key to its applicability to spatially extended systems. 

Keywords. Bayesian inference, Monte Carlo method, sequential data assimilation, ensemble Kalman filter, 
particle filter, localization, dynamical systems. 
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1 Introduction 

This paper is concerned with algorithms for assimilating observational data into computational models using 
ensemble prediction techniques and sequential Bayesian model updates. The focus will be on state estimation 
problems which give rise to the classic filtering problem of stochastic processes El- 

There are two popular classes of ensemble based filtering techniques; namely particle filters (PFs) [5] and 
ensemble Kalman filters (EnKFs) [10]. While PFs lead to consistent approximations of the underlying filtering 
problem [5], they suffer from the curse of dimensionality |1|. EnKFs, on the other hand, are applicable to 
spatially extended systems when combined with inflation and localization nni; but they lead to estimates which 
are inconsistent with the Bayesian approach to filtering unless all involved probability density functions (PDFs) 
are Gaussian. There have been several recent attempts to enhance the performance of PFs by either using 
better proposal densities (see [31j for a review) or by implementing localisation into PFs |2ll[6|. At the same 
time, several efforts have been made to extend EnKFs beyond their Gaussian prior assumption. Most of these 
attempts are built around Gaussian mixture approximations and, hence, can be seen as efforts to bridge EnKFs 
with PFs (see, for example, UHlIMlIiaill])- Other recent attempts to extend EnKFs beyond Gausianity include 

umanii. 

In this paper, we combine the general linear ensemble transform filter (LETF) framework, as laid out in detail 
in m, with the EnKF-PF bridging approach suggested in [11]. As in m, our hybrid ensemble filter approach 
bridges EnKFs with PFs through a parameter a S [0,1], which can be chosen adaptively. Key innovations of 
our approach, as compared to the hybrid method of include the use of ^-localisation [Tl[17], the fact that 
either a PF or an EnKF can be applied first, that any implementation of an EnKF can be used (including [T]), 
and that observation operators need not to be linear. 

The properties of our hybrid ensemble filter will be demonstrated numerically for the Lorenz-63 model 
[T5] . the Lorenz-96 model [MUD], and a coupled Lorenz-96 wave equation model The numerical results 
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demonstrate that our hybrid approach leads to reduced root mean square (RMS) errors compared to an EnKF 
even for ensemble sizes at which a standard PF fails to track the reference solution. 

We now summarize the basic data assimilation problem which will be considered throughout this paper. See 
nanz] for further background material on filtering and general data assimilation algorithms. The basic task 
of sequential filtering (or sequential data assimilation) is to estimate a reference trajectory Zref{t) G of an 
evolution equation, such as 

|=/(M). (1) 

from partial and noisy observations 

Vohsij'k) ~ ^(■^ref (^fc); ^) “f A; = 1, 2, ... , (2) 

at discrete times tk > 0. Here h : x M>o —>■ is the observation operator and G fc > 1, denote 

independent realizations of a Gaussian random variable with mean zero and covariance matrix R G 
We are also given the PDF, 7rz(z,0), for the initial conditions z(0) of 0 at time zero. In order to simplify 
the subsequent discussion, we will assume that the observation operator h is linear and time independent, i.e., 
h{z, t) = Hz. 

FnKFs and PFs both aim at providing approximations to the conditional distributions 

TTziz, t|j/obs(Al), • • • , Vohsitk)) 

for t > t/c in the form of an ensemble Zi{t), i = 1,... ,M. The best estimate for Zref(A) is provided by the 
ensemble mean, 

1 “ 

^est{i) = ^ ^ (i-) ■ (3) 

^ i=l 

While Zest (A) from a PF converges to the true expectation value 

z{t)= znz{z,t\yohs{ti),-■ ■ ,yohs{tk))dz (4) 

as M —>■ oo, the same does not hold true for FnKFs, in general. However, it is often found in practice that 
FnKFs lead to smaller time-averaged RMS errors, 

RMSF=^^ y^||Zest(A) - Zref(i)Pdt, (5) 

for small or moderate ensemble sizes, M. See [8] for more details. 

Any ensemble-based filter algorithm alternates between prediction steps, where an ensemble is propagated 
under the dynamics 0, and analysis steps, where a forecast ensemble zf = Zi{tk) is transformed into an analysis 
ensemble zf based on the available observations at time t = tk- The analysis ensemble then provides the initial 
conditions for the next prediction step over the time interval t G It is the implementation of the 

analysis step that provides the crucial difference between the various filter algorithms. Section lays out a 
general framework for formulating analysis steps and discusses the ensemble square root filter (FSRF) [30] and 
the ensemble transform particle filter (FTPF) |5] as two specific examples for FnKFs and PFs, respectively, 

which will be used throughout this paper. The novel hybrid filter is based on a combination of the FSRF and 

the FTPF and is introduced in detail in Section |2.3[ Its generalization to spatially extended systems using 
localization is formulated in Section The complete hybrid algorithm is summarized in Section [^ while the 
numerical results can be found in Section jsj 


2 Summary of the LETF framework 

Given a forecast ensemble z^, i = I,..., M, at an observation time tfc, we consider general linear transformations 


M 


= 


( 6 ) 


2 



into an analysis ensemble 2 ;^", j = 1,..., M. The coefficients dij satisfy 

M 

(7) 

i=l 

and depend on the forecast ensemble and the observed j/obs = 2/obs(ifc)- We demonstrate in the next two 
paragraphs how the ESRF and the ETPF fit into the LETF framework. See [57] for further details and LETF 
formulations of other EnKFs in particular. 


2.1 ESRF formulation 

We introduce the forecast ensemble mean 




1 ^ 


2=1 


the Nz X M matrix of ensemble anomalies 

klf=[(4-Z-f) (4-Z-f) ... 

the symmetric M x M matrix 


s= U + --J—{haYr~^ha^ 

' M -1^ ' 


- 1/2 


and the weight vector w G with entries 

1 1 


Wi = - 

* M M- 


— e/ S\HA^YR-\H-z^ - yobs). 

Here Ci denotes the zth basis vector in Also note that the entries Sij of the matrix S satisfy 


MM M 

T ^ X! ^ ^ y] Wi = 1. 

^=1 2=1 2=1 

Following the presentation in m, the analysis ensemble of an ESRF is now given by 


M 


M 


= E = E 4d: 


with filter coefficients 


iy=rfy ({4} 


( 8 ) 

(9) 

( 10 ) 

(11) 

( 12 ) 


M 

2=1 

(13) 


(14) 


which depend on the forecast ensemble and the observation, T/obs- We note that the analysis mean satisfies 

M 


= E 


Z^,Wi. 


(15) 


i=l 


In case of nonlinear observation operators, one would replace HA^ by the Ny x M matrix 

[{h{z{,tk)-hY {h{z\Ak)-h^) ■■■ {h{zYAk)-hY] ( 16 ) 

and Hz^ by 

1 ^ 

(17) 


i=l 


in the formulas above. We finally mention that the computational complexity of the ESRF scales like 0{M^) 
in the ensemble size M. 
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2.2 ETPF formulation 


The ETPF provides an implementation of PFs which replaces the resampling step of a sequential importance 
resampling (SIR) PF by a linear transformation of the form ([^ [211 [3 HZ]. As for a SIR PF, we define 
importance weights 

exp - yohs)^R~^{Hzl - j/obs)) 


Wi = 


Y.j=i exp {-\{Hz^. - VohsVR~^{Hz^. - 2 /obs)) 


The required transformation ®> is then found by solving the linear transport problem 


M 


T* = arg min t. 


where 11 denotes the set oi M x M matrices T with non-negative entries tij subject to 

M M 

tij = 1, and tij = WiM. 

i=i j=i 

Finally, the ETPF defines the analysis ensemble via 


(18) 


(19) 


( 20 ) 


M 


Rj 


with filter coefficients 
and analysis ensemble mean 


4 = 

dW = d^4({4}^yohs) ■■= t*i 


M MM M 

= M E y = M E E ^ 4 , = E 

j — l 3 — ^ i—1 i—1 


( 21 ) 

( 22 ) 

(23) 


The consistency of a single ETPF step with Bayes’ formula has been demonstrated in [21] ■ However, the ETPF 
can fail to track the reference solution, Zref(t), if the effective sample size 


MeB = 


EIVl 2 


(24) 


becomes too small. In the extreme case of M^b = !> all analysis ensemble members, Zj, become equal to zj,, 
where z^, is the forecast ensemble member with importance weight Wi* = 1. This limitation of the ETPF 
provides the main motivation for our hybrid ensemble transform particle filter. 

The linear transport problem (19) can be solved by a linear programming algorithm [22] . In our numerical 
experiments we used the FastEMD algorithm of [22] j which has a computational complexity of 0{M^ logM). 
We also note that one dimensional transport problems can be solved efficiently by sorting the ensemble members 

HD- 

Nonlinear observation operators can be accommodated by the ETPF by simply replacing Hzl by h{zl,tk) 
in the computation of the importance weights (18). 

2.3 The hybrid ensemble transform particle filter 

Following HD, we split the likelihood function into a product of two factors, i.e., 

7ry(j/obs|2) oc exp (^-‘^{Hz - yohs4R~^{Hz - yobs)) x exp ^ ^ {Hz - yohs4R~^{Hz - yobs)^ (25) 
for a G [0, Ij. The second factor is treated by the ESRF giving rise to 

- 1/2 




(26) 
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( 27 ) 


weight vector w{a) G with entries 


and filter coefficients 




d^^{oiAz!},yohs) ■■= m{a) - +Sij (a)- 


The first factor is treated by the ETPF with importance weights 

exp (-f (i7zf - - 2 /obs)) 

'^2(0^) = -T?-^-=- - - 

exp - yohsVR~^{Hz] - yobs)) 


and the linear transport problem 


M 


T*{a) = argmin t 




i j = i 


where 11 now denotes the set oi M x M matrices T with tij > 0 subject to 

M M 

= 1, and tjj = Wi{a) M. 

i=l j=l 

The coefficients of the associated ETPF are then given by 

dWi(^,{zl},yOhs) ■■= t*Aa). 


(28) 

(29) 

(30) 

(31) 

(32) 


When implementing the hybrid filter, one can decide between two options; either one applies the ETPF or 
the ESRF hrst. The two resulting implementations can be summarized as follows: 

(A) ETPF-ESRF 


(B) ESRF-ETPF 


M 

i—l 

“y ■ 

= dfj^{a,{zl},yohs), 

M 


■.= d^^{a,{zn,yohs). 

M 

Y.-'xW. 

i—l 

d^^ 

■= d^^{a,{zj},y Ohs), 

M 

2=1 

dA^ 

:=4^(a,{zf},yob,). 


(33) 

(34) 


(35) 

(36) 


Note that (341 uses the updated ensemble from (33) to determined the coefficients d^^. The same dependency 
applies to the ESRF-ETPF. 

The bridging parameter a can either be set to a constant value or selected adaptively El such that the 
ratio between the effective sample size, 


Mes{a) := 


1 


(37) 


in the ETPF part of the hybrid filter and the ensemble size, M, is equal to a fixed reference value 9 G [0,1], i.e. 

(38) 


Note that 6 = 1 implies a = 0 and M^ef{oi) is monotonically decreasing in a. If M^^Aa) > 9M for all a G [0,1], 
then we set a = 1. 
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2.4 Ensemble inflation and particle rejnvenation 

Filter implementations with small or moderate ensemble sizes often lead to filter divergence, i.e., the filter is 
no longer able to track the reference solution. Filter divergence can often be traced back to analysis ensembles 
becoming underdispersive, i.e., the analysis ensemble spread does not properly represent Bayesian posterior 
uncertainties. In order to prevent ensembles from becoming underdispersive, one can apply multiplicative 
ensemble inflation to the forecast ensemble prior to an analysis step and/or one can apply particle rejuvenation 
to the analysis ensemble (see, for example, [101 HZ]). We used particle rejuvenation in the form of 


M 

2=1 






y/M-I 


(39) 


in our numerical experiments. Here /3 > 0 denotes the rejuvenation parameter and the are i.i.d. Gaussian 
random variables with mean zero and variance one. We also enforce 

M 

E^b = 0 (40) 


in order to preserve the ensemble mean, under rejuvenation. 

Note that the particle rejuvenation step is the only instance in the complete filter algorithm where random¬ 
ness enters. This randomness provides a safeguard against the creation of identically analysis ensemble members 
under an ETPF, as discussed towards the end of Section [2^ 


3 Spatially extended systems and localization 

At each data assimilation step, we have to now deal with ensembles of spatial fields z(x) G where x G 
denotes the spatial variable of dimension d. Those fields are provided on a computational grid with grid points 
Xk G k = 1,... ,K. Spatially discrete fields are formally collected into a state vector 

z = (zixi)'^ z{x 2 )'^ ■■■ z{xk)'^)'^ (41) 

of dimension Nz = N x K. A forecast ensemble 

4 = (zUx.f zlix^f ■ ■ ■ zlixKff G (42) 

z = 1 ,..., M, is then transformed at each grid point, according to 

M 

^ ^ (x/e); }':yohs}i ('^^) 

2=1 

which gives rise to the analysis ensemble 

Zj = izjixi)'^ Zj{x2)'^ ■ ■ ■ Zjixx)'^)'^ G (44) 

j = 1,... ,M. It is crucial that these local updates are implemented in a manner that maintains the spatial 
regularity of the forecast states. This prevents the use of localization in a standard SIR PF. 

We assume that observations are taken at a discrete subset of grid points q = 1,..., Q, and satisfy the 
observation model 

Vq = Uzilq) + r]q 

with independent measurement errors rjg ~ N(0,r) and observation operator H : -G We assume here, 

for simplicity, Ny = 1 and r > 0 denotes the variance of the measurement error. We introduce the more compact 
notation 

Vohs = Hz+ 11 (45) 

with Uohs = (yi,... ,j/Qy G M‘5, r] = ( 771 ,..., Tyg)’’’ ~ N(0, i?), R = rl, and H : -G M'S appropriately 

defined. If each grid point is observed then we simply have Q = K, fg = Xg with q = 1,... ,K. If every second 
grid point is observed then Q = K/2 and y, = X 2 q with q = 1,... ,K/2 etc. 
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We now extend the concept of i?-localization [HISIIIT] to the hybrid filters introduced in Section [2^ First 
we define a localized inverse measurement error covariance matrix R{xk)~^ & with diagonal entries 






(46) 


q = for given localization radius i?ioc > 0 and compactly supported tempering function p{t) with 

p(0) = 1 and p{t) = 0 for t > 2. We used 


0 


1 - + 1^3 + 1^4 _ 1^5 foj. ^ 

p{t) = + 4 — 5t + for 1 < t < 2, 


(47) 


otherwise , 


in our numerical experiments [12) . The two hybrid algorithms from Section 2.3 are now adjusted as follows. 
The ESRF is replaced by the LETKF [13], which leads to 


S{xk, ak) = U+ \:^{HA^fRixk)-^HA^ 


- 1/2 


Wiixk,ak)-^ M-1 


M-1 

^^eJS{xk,akf‘{HA^fR{xk)~^{Hz^-yohs), 


and 


^ij }; 2/obs) • 


(48) 

(49) 

(50) 


The ETPF filter contribution is modified at each grid point Xk as proposed in [3[17|. Specifically, we define 
localized weights 

exp - yohs)'^R{xk)~^{Hzl - yohs)) 


exp (-^(i 7 zf - yohsVRixk)~^{Hz^^ - yohs)) 

introduce the linear transport problem 


(51) 


M 

T* = arg min tij\\zl{xk) - z)i{xk)\\^ (52) 

ij'=l 

at each grid point, Xk, where 11 denotes the set oi M x M matrices T with tij > 0 subject to 

M M 

kj = 1, and kj = Wi{xk, Ofc) M, (53) 

i=i 


and define 


^ij (^fc ; Cl/e, {2:^ }, ^obs) — ^ 


ij ■ 


(54) 


We note that the distance term in (521 can be replaced by other distance measures involving, for example, grid 


point values in the vicinity of Xk- See laisT] for more details. The cost functional ( |5^ has the advantage that 
the associated transport problem is easy to solve in case z{xk) G R, i.e. = 1. 

Finally, the desired localized hybrid filter schemes can be implemented in the following two variants: 

(A) ETPF-LETKF 


M 


z)^{xk) =^zl{xk)dff (xk), dff {xk) := dff {xk,ak,{zf},yohs), over all grid points, (55) 

M 

Zj{xk) = X! ^ti^k)d^^{xk), df^{xk) ■■= d!ff{xk, ak, {z^}, 2/obs), over all grid points. (56) 


7 











(B) LETKF-ETPF 


M 

Zj{xk) = ^ zl{xk)d^^{xk), df^{xk) ■= df^{xk, au, {zf}, j/obs), over all grid points, (57) 

M 

Zj{xk) = ^ Zi{xk)d^j^{xk), d^f{xk) ■■= df/(a:fc, ak, { 2 ?}, yobs), over all grid points. (58) 

i=l 


The bridging parameters, ak, can again be either set to a constant value, a, for all grid points or determined 
adaptively. In the latter case ak is determined by 


M^e{{xkak) 

M 


— 0 , ■ — 


Yl,fiiWY{xk,aky 


(59) 


for given effective ensemble size ratio, 9 G [0,1]. If MYef{xk, ak) > 9M for all ak G [0,1], then we set a*, = 1. 
Finally, particle rejuvenation in the form of (|39|) is performed at the end of each filter step. 


4 Summary of hybrid filter algorithm 


We summarize the numerical implementation of the proposed hybrid filter as it has been used for the numerical 
experiments below. We explicitly describe the filter implementation with localization. An implementation 
without localization follows from straightforward modifications. The free parameters of the data assimilation 
steps are the ensemble size, M > 1, the localization radius, i?ioc > 0, in (461, the bridging parameters, ak G [0,1], 
which are either fixed to a value a or determined adaptively through the effective sample size ratio, 9 G [0,1], 
in (59). Any choice of these parameters together with a forecast ensemble Z = 1,..., M, and an observation 
VohAt^k) determine the ETPF-LETKF update step (55l-(56) uniquely. The linear transport problems, defined 
by (52 )-(|5^ are solved using the Matlab implementation of FastEMD [53]. 


The update step (55)-(56) is followed by particle rejuvenation (39) with parameter /3 > 0. Note that 
rejuvenation is applied to the whole spatially dependent analysis fields with one and the same set of random 
numbers ^ij, i.e. 


zj(a:fc) 


zKxk) 




M 

Y^{zl{xk) - z\xk))--^:d^ 

i—1 ^ 


(60) 


for all grid points Xk- 

A complete implementation also requires a time-stepping method for the differential equation 0- We used 
the implicit midpoint rule with step-size denoted by At > 0. Observations ([^ are generated from a reference 
trajectory, {z"gf}„>o, which is computed using the same time-stepping method in ([^ and a randomly selected 
initial condition 2 )),^ ^ nz{z,Q)- The initial ensemble, z^, i = 1,...,M, is also drawn from the initial PDF, 
ttz{z,Q). 


5 Numerical experiments 

We now conduct a series of numerical experiments in order to assess the relative merits of the hybrid approach 
compared to a purely EnKF or PF based data assimilation approach. The ensemble sizes, M, are chosen such 
that standard PF implementations cannot compete with an EnKF. Such a data assimilation scenario is typical 
in complex application areas such as numerical weather prediction m- The quality of a filter is assessed by 
the following time-discretised version of ([^: 

1 ^ n- 

RMSE = - ^ J - ZYetitk)P. (61) 

Here K denotes the total number of assimilation steps. 
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Figure 1: Hybrid filter for Lorenz-63 system with ensemble sizes M = 15,... ,35: left panel: RMS errors for 
fixed bridging parameter; right panel: RMS errors for adaptively determined bridging parameter. Note that 
a = 0 and 9 = 1 correspond to a pure ESRF filter. 


5.1 Lorenz-63 


The first numerical tests of our hybrid filter are performed on the chaotic Lorenz-63 system m with its standard 
parameter setting a = 10, p = 28, and /3 = 8/3. We only observe the first component of the three dimensional 
system in observation intervals of Afobs = 0.12, i.e., tk = 0.12A: in (611, and with observation error variance 


R = 8. A total oi K = 100,000 assimilation steps are performed and = 3 in (61). The differential equations 
are solved numerically with a step-size of At = 0.01. 

The ETPF-ESRF filter is implemented for ensemble sizes varying between M = 15 and M = 35. Particle 
rejuvenation is applied with /? = 0.2. Simulations with /3 = 0.15 and /? = 0.25 gave similar results to those 
reported here. No localization is necessary for this low dimensional system. Our numerical results for fixed and 
adaptively determined bridging parameter, a, can be found in Figure]^ It can be seen that a fixed bridging 
parameter of a = 0.2 leads already to substantial improvements for the smallest ensemble size of M = 15. 
Further improvements are achieved for larger ensemble sizes at larger bridging parameters. In contrast to these 
very encouraging results for fixed bridging parameters, the adaptive implementation of the hybrid filter does 
not work quite as effectively for small ensemble sizes of M = 15 and M = 20. This finding indicates that one 
might have to search for other selection criteria. 

We also implemented the Ensemble Kalman Particle Filter of m for comparison. The Ensemble Kalman 
Particle Filter bridges the EnKF with perturbed observations [101 ^ Here we replaced the SIR PF 

with the ETPF in order to facilitate a fair comparison. The numerical results are displayed in Figurej^ We find 
that the Ensemble Kalman Particle Filter leads to larger RMS errors for all values of the bridging parameter, 
a. We speculate that this is due to the use of the EnKF with perturbed observations and the fact that the 
EnKF is applied first. This hypothesis is supported by the numerical results from the ESRF-ETPF, which also 
does not lead to significant improvements for a > 0. See Figure 

Note that we run the hybrid filters at ensemble sizes, M, for which EnKFs outperform PFs. The more 
interesting it is that the ETPF-ESRF leads to substantial improvements for values of the bridging parameter 
in the range [0.2, 0.4] for ensemble sizes between M = 15 and M = 35. Figure shows that the optimal a 
increases with increasing ensemble sizes, M. Asymptotically, one expects that a —>■ 1 and 0 —^ 0 as M —>■ oo. 
Indeed, we found numerically that a = 0 is optimal for ensemble sizes smaller than M = 10 and that a = 1 for 
ensemble sizes larger than M = 100. See also the convergence study in Section [5^ 


5.2 Lorenz-96 

We next test our hybrid approach on the spatially extended Lorenz-96 model 

±1 = {xi+i-xi-2)xi-i-xi+8, Z = l,...,40 (62) 
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Ensemble Kalman Particle Filter hybrid ESRF-ETPF 



Figure 2: Left panel: RMS errors for Ensemble Kalman Particle Filter applied to the Lorenz-63 system with 
ensemble sizes M = 15,... ,35 and for fixed bridging parameter; right panel: RMS errors for the ESRF-ETPF 
using the same setting and a in the range [0,0.5]. 


[ini [20]. We observe every other grid point in time intervals of Atobs = 0.11 with independent Gaussian 
measurement errors of variance equal to eight. A total oi K = 50, 000 assimilation steps is performed. The 
hybrid ETPF-LETKF filter is implemented with localization radius i?ioc = 4, particle rejuvenation /3 = 0.2 and 
ensemble sizes M = 20, 25, 30. The equations of motion are integrated numerically with the implicit midpoint 
rule and At = Atobs/22. The resulting time-averaged RMSE can be found in Figurej^ The overall findings are 
similar to the Lorenz-63 model except that the adaptive approach now performs comparably (or even better) to 
an optimally tuned fixed parameter implementation and that the overall improvements in terms of RMS errors 
are less significant. The later fact is likely due to the weaker nonlinearity of the Lorenz-96 model. We also 
implemented the hybrid filter for localization radii i?ioc G {3, 5, 6,7,8} and found that a > 0 always leads to 
improvements in the RMS error compared to a pure LETKF implementation under the simulation conditions 
stated above. We finally implemented the Lorenz-96 in the more traditional setting of each grid point being 
observed in intervals of Atobs = 0.05 with measurement error variance r = 1. Our numerical findings indicate 
that the achievable improvements of our hybrid filter in comparison to the LETKF are rather minor and require 
ensemble sizes M > 30. We believe that these hndings reflect the fact that frequent observations of the complete 
state space lead to nearly Gaussian forecast distributions. 


5.3 Lorenz-96 coupled to a wave equation 

Physical processes in the geosciences often act on different time-scales. To simulate that effect we use a model, 
that consists of the advective Lorenz-96 system (62) and relatively fast moving waves, which are modeled by 
the discrete wave equation 


e^hi ——hi + c^{hi+i — 2hi + hi-i) — 2e^^hi, Z — 1,...,40. (63) 

The positive parameter e ^ 1 characterises the ratio between the advective velocity in x and the wave speed 
of h. The parameter c > 0 determines the wave dispersion over the computational grid and 7 > 0 is a damping 
parameter. The coupled model with coupling-strength 6 looks as follows [5|: 

xi = {l- S){xi+i - xi_ 2 )xi-i + S{xi_ihi+i - xi_ 2 hi-i) -xi + 8, (64) 

e^hi = -hi + c^{hi+i - 2hi + hi^i) + xi - 2e^jhi, (65) 

1 = 1,..., 40, and periodic boundary conditions are applied. There exists a balance relation 

xi = hi- c^{hi+i - 2hi + h;_i) ( 66 ) 


between the x and h variables. If it is fulfilled initially, solutions stay in approximate balance for long time 
intervals with the deviations from exact balance proportional to e. For our experiments, we used this model 
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Figure 3: Hybrid filter for Lorenz-96 system with ensemble sizes M = 20, 25, 30: left panel: RMS errors for 
fixed bridging parameter; right panel: RMS errors for adaptively determined bridging parameters 


with a coupling strength of <5 = 0.1, e = 0.0025, a damping of 7 = 0.1, and c = 1/2. The time-stepping method 
of [S] is used with a step-size of At = 0.002. 

The initial ensemble in {a;;} is created by sampling from a Gaussian distribution with mean zero and 
covariance 0.01/. The /-values of the initial ensemble are created by calculating them from the x-values using 
the balance relation ( 661 . The /-values are set equal to zero. We observe xi at every other grid-point at time 
intervals of /Siobs = 0.15 with measurement-error variance R = 8 / 20 . We applied the hybrid ETPF-LETKF 
filter to our model with a localization radius of R\oc = 4 and a particle-rejuvenation parameter of /3 = 0.2. We 
let the experiment run for K = 50,000 assimilation cycles and a spin-up phase of 1000 cycles, in which a is set 
to zero and the /-values are calculated to be in exact balance to the x-values after every assimilation step. We 
conduct experiments for ensemble sizes of M = 20,25,30, fixed bridging parameters a = 0,0.1, ...,0.6 and also 
for adaptive Ofc’s using 9 = 0.75,0.8,..., 1. 

The resulting time-averaged RMS errors are displayed in Figure]^ One can see that the minimal RMS error 
is obtained at a = 0.3 when M = 20, at a = 0.4 when M = 25 and at a = 0.5 when M = 30. The results 
for the adaptive approach show, that the RMS error is minimal for values of / G [0.90,0.95] depending on the 
ensemble size, M, and that the adaptive approach is less efficient than an optimal choice of a. This is likely due 
to the fact that an adaptive choice of the bridging parameter in each grid point, Xk, together with localization 
leads to less balanced fields. 

We finally mention that the damping in (65) is necessary to suppress the accumulation of imbalances gener¬ 
ated in each assimilation step. See for a detailed discussion of this phenomenon and a potential cure in the 
form of a mollified filter approach. 


5.4 Convergence stndy for a single Bayesian inference step 

After conducting a series of numerical experiments for which the ensemble sizes were relatively small, we finally 
demonstrate that the proposed hybrid filter indeed bridges EnKFs and PFs as the ensemble size varies from very 
small to very large. In order to keep the computations feasible, we only perform a singe assimilation step with 
given bimodal prior (see Figure]^ and the ETPF-ESRF implemented for ensemble sizes M G [2,4, 8 ,..., 1024]. 
Each experiment is repeated 100,000 times and the averaged RMS errors for the estimated posterior mean are 
computed in terms of M and a. The optimal value of a and its RMS error are displayed in Figure[^as a function 
of the ensemble size, M. As expected, the optimal value of a is close to zero for small M and approaches one 
as M increases. 
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Figure 4: Hybrid filter for Lorenz-96 coupled to a wave equation with ensemble sizes M = 20,25, 30: left 
panel: RMS errors for fixed bridging parameter; right panel: RMS errors for adaptively determined bridging 
parameters 



ETPF-ESRF 




Figure 5: Single data assimilation step with ETPF-ESRF; left panel: prior and posterior distributions; right 
panel: optimal bridging parameter and resulting RMS errors as a function of the ensemble size. 
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6 Conclusions 


Further numerical experiments show that the ETPF-ESRF and ETPF-LETKF, respectively, outperform the 
ESRF-ETPF and LETKF-ETPF, respectively, on all numerical examples of Section]^ See [7j for more details. 
We believe that this is due to the fact that the true forecast/prior distributions in our examples are more 
non-Gaussian than the analysis/posterior distributions. This is the case, in particular, for the Lorenz 63 model. 
Under those circumstances it seems preferable to apply a PF first followed by an EnKF. We also found that 
the adaptive hybrid method works well for the Lorenz-96 system, while a fixed choice of a is preferable for the 
Lorenz-63 model. It is feasible that other adaptation criteria for the choice of a could lead to better results. 
For example, one could select a such that the Kullback-Leibler divergence between the uniform measure, 1/M, 
and the importance weights (29) remains below a given threshold. 

Our numerical results also indicate that our hybrid method is suitable for multi-scale systems. However, 
our hybrid method does not address the problem of generating unbalanced fields through localization (see, for 
example, 0 ). This issue will be addressed in a separate publication using the idea of mollification, as proposed 
in [^. We finally mention that our hybrid methods can be combined with modihed proposal densities for PFs 
(see [3T] for a recent review of the subject). 
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